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i -^h ■ We describe a mathematical formalism and numerical algorithms for identifying 

and tracking slowly mixing objects in nonautonomous dynamical systems. In the au- 
tonomous setting, such objects are variously known as almost-invariant sets, metastable 
sets, persistent patterns, or strange eigenmodes, and have proved to be important in 
£Nj . a variety of applications. In this current work, we explain how to extend existing au- 

tonomous approaches to the nonautonomous setting. We call the new time-dependent 
slowly mixing objects coherent sets as they represent regions of phase space that dis- 
perse very slowly and remain coherent. The new methods are illustrated via detailed 
examples in both discrete and continuous time. 

1 Introduction 



The study of transport and mixing in dynamical systems has received considerable attention 
in the last two decades; see e.g. [311 ES H HZ] for discussions of transport phenomena. In 
particular, the detection of very slowly mixing objects, known variously as almost-invariant 
sets, metastable sets, persistent patterns, or strange eigenmodes, has found wide application 
in fields such as fluid dynamics [3H 1301 ES], ocean dynamics (2TJ ?], astrodynamics [I], and 
molecular dynamics (HUH]. A shortcoming of this prior work, based around eigenfunctions of 
Perron-Frobenius operators (or transfer operators, or evolution operators) is the restriction 
to autonomous systems or periodically forced systems. In this work, we extend the notions of 
almost-invariant sets, metastable sets, persistent patterns, and strange eigenmodes to time- 
dependent Lagrangian coherent sets. These coherent sets form a time parameterised family 
of sets that approximately follow the flow and disperse very slowly; in other words they stay 
coherent. Coherent sets are the natural nonautonomous analogue to almost-invariant sets. 

The standard dynamical systems model of transport assumes that the motion of passive 
particles are completely determined by either an autonomous or a time-dependent vector 
field. Traditional approaches to understanding transport are based upon the determination 
of the location of geometric objects such as invariant manifolds. In the autonomous setting, 
an invariant manifold of one dimension less than the ambient space will form an impenetrable 
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transport barrier that locally partitions the ambient space. In the periodically-forced setting, 
primarily in two-dimensional flows, it has been shown that slow mixing in the neighbourhood 
of invariant manifolds is sometimes controlled by "lobe dynamics" (361 E3 HE] . In the truly 
non-autonomous, or aperiodically forced setting, finite-time hyperbolic material lines [21] 
and surfaces [25] have been proposed as generalisations of invariant manifolds that form 
barriers to mixing. These material lines and surfaces are known as Lagrangian coherent 
structures; see also [42 J for an alternative definition. The geometric approach can often be 
used to find co-dimension 1 sets (coherent structures) that form boundaries of coherent sets. 

An alternative to the geometric approach is the ergodic theoretic approach, which at- 
tempts to locate almost-invariant sets (or metastable sets) directly, rather than inferring 
their location indirectly from their boundaries. The basic tool is the Perron-Frobenius op- 
erator (or transfer operator). Real eigenvalues of this operator close to 1 correspond to 
eigenmodes that decay at slow (exponential) rates. Almost-invariant sets are heuristically 
determined from the corresponding eigenfunctions / as sets of the form {/ > c} or {/ < c} 
for thresholds cGK. Such an approach arose in the context of smooth autonomous maps 
and flows on subsets of M d E] about a decade ago. Further theoretical and computational 
extensions have since been constructed [T7J [151 HB]- A parallel series of work specific to 
time-symmetric Markov processes and applied to identifying molecular conformations was 
developed in [IHl El [9] and surveyed in [4T] . 

There have been some recent studies of the connections between slow mixing in period- 
ically driven fluid flow and eigenfunctions of Perron-Frobenius operators. Liu and Haller 
[30] observe via simulation a transient "strange eigenmode" as predicted by classical Flo- 
quet theory. Pikovsky and Popovych [33] [35] numerically integrated an advection-diffusion 
equation to simulate the evolution of a passive scalar, observing that it is the sub-dominant 
eigenfunction of the Perron-Frobenius operator that describes the most persistent deviation 
from the unique steady state. 

The Perron-Frobenius operator based approach has been successful in a variety of appli- 
cation areas, however, as the key mathematical object is an ezgenfunction, there is no simple 
extension of the method to systems that have nonperiodic time dependence^ 

Indeed, Liu and Haller [3U] state that: 

"...strange eigenmodes may also be viewed as eigenfunctions of an appropri- 
ate Frobenius-Perron operator... This fresh approach offers an alternative view 
on scalar mixing, but leaves the questions of completeness and general time- 
dependence open." 

It is this question of general time-dependence that we address in the current work. We 
extend a standard formalism for random dynamical systems to the level of Perron-Frobenius 
operators to create a Perron-Frobenius operator framework for general time-dependence. 
We also state an accompanying numerical algorithm, and demonstrate its effectiveness in 
identifying strange eigenmodes and coherent sets. 

1 A relevant analogy to see this is the following. Consider repeated application of a single matrix A. 
The eigenvectors of A provide information on directions of exponential growth/decay specified by the cor- 
responding eigenvalues. Similarly, the eigenvectors of a product of matrices Ak ■ ■ ■ A2A1 describe directions 
of exponential growth/decay, specified by the eigenvalues of the product, under repeated application of 
this matrix product. However, the directions of exponential growth/decay under a non-repeating product 
• ■ ■ Ak ■ ■ ■ A 2 Ai cannot be in general be found as eigenvectors of some matrix. 
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An outline of the paper is as follows. In Section 2 we formalise the notions of nonau- 
tonomous systems in both discrete and continuous time. In Section 3 we describe a Galerkin 
projection method that we will use to produce finite matrix representations of Perron- 
Frobenius operators. In Section 4 we define the critical constructions for the nonautonomous 
setting. We show that the nonautonomous analogues of strange eigenmodes are described 
by the "Oseledets subspaces" or "Lyapunov vectors" corresponding to compositions of the 
projected Perron-Frobenius operators. In Section 5 we describe in detail a numerical al- 
gorithm to practically compute these slowly decaying modes, and demonstrate that in the 
continuous time setting, these modes vary continuously in time. Our numerical approach 
is illustrated firstly in the discrete time setting with an aperiodic composition of interval 
maps, and secondly in the continuous time setting with an aperiodically forced flow on a 
cylinder. Section 6 provides some further background on almost-invariant sets and coherent 
sets and Section 7 describes a new heuristic to extract coherent sets from slowly decaying 
modes in the nonautonomous setting. This heuristic is then illustrated using the examples 
from Section 5. 



2 Nonautonomous Dynamical Systems 

We will treat time dependent dynamical systems on a smooth compact <i-dimensional mani- 
fold M C K D , D > d in both discrete and continuous time. In order to keep track of "time" 
we use a probability space (O, "H, P), with the passing of time controlled by an ergodic auto- 
morphism 9 : Q O preserving P (ie. P = Po#~' for all t > 0). We require this somewhat more 
complicated description of time for technical reasons: to run the ergodic-theoretic arguments 
in Theorem [1] that guarantee the existence of the nonautonomous analogues of strange eigen- 
modes. The requirement that P be an ergodic probability measure rules out obvious choices 
for Q and 9: (i) in discrete time, Q = Z and 9 s (t) — t + s, and (ii) in continuous time, Q = R 
and 9 s (t) = t + s. In both (i) and (ii), there is no ergodic probability measure on Q preserved 
by 9. In the next two sections, we will introduce suitable examples of Q and 9 and describe 
the nonautonomous systems they generate. 

2.1 Discrete time — Maps 

In the discrete time setting, we will think of Q C (Z) z , and 9 as a left shift a on Q defined 
by (<ru)i = where u = (. . . , cj , cji, . . .) G Q. We assume that a is ergodic with 

respect to P. Let T = {T^,,}^^ be a collection of (possibly non-invertible) piecewise 
differentiable maps on a compact manifold M. For brevity, we will sometimes write T w in 
place of T^. We will define a nonautonomous dynamical system by map compositions of 
the form T a k-\^ o • ■ ■ o T auJ o T u . Define 

{Tak-i^o ■ ■ ■ oT^oT.fi), k > 0; 
Id, k = 0; 

T~} k o-oTl oT:\ (x), k<0. 

For k > (resp. k < 0), $>(k,uj,x) represents the forward time (resp. backward time) fc-fold 
application of the nonautonomous dynamics to the point x initialised at "time" u. Whenever 
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G<t>Q 



Figure 1: Graph of the sequence space fi. 



T u is non-invertible, T w will represent the finite set of all preimages of x. We call $ a 
map cocycle. 

Definition 1. Endow M with the Borel a-algebra and let [i be a probability measure on 
M. We call /i an invariant measure if /i o $(— 1, ■) = fi for all Co> G fi. 

This definition of an invariant measure is stricter than is usual for random or nonau- 
tonomous dynamical systems (e.g. (2j Definition 1.4.1]). More generally, one may allow 
sample measures jj, = n u and insist that ° 1, •) = fi w for all u E Q. 

Example 1 (Aperiodic map cocycle). We construct a map cocycle $ by the composition 
of maps Ti from a collection T according to sequences of indices u 6 H. The collection 
T := {T 1 ,T 2 ,T 3 ,T 4 } consists of expanding maps of the circle S 1 , which we think of as [0, 1] 
with endpoints identified. The sequence space Q C {1,2,3,4} Z is given by 



n = {w6{l 1 2,3, 4} Z : e Z, M WiU)i+1 = 1 } 



with adjacency matrix 



M 



( 1 

1 



1 



1 / 



Elements of Q correspond to bi-infinite paths in the graph Figure 1. The shift o : Q — > Q 
is a subshift of finite type. A Borel a -algebra % is generated by the length-one cylinder sets 
Ci = {u : uq = i} , i = I, . . . ,4, and by giving equal measure to these four cylinder sets, we 
generate a shift-invariant probability measure P. 

The maps ofT are defined in terms of a continuous piecewise-linear map H a : S 1 — > S 1 , 
which has almost-invariant sets (see Definition^ Section 6) [0,0.5] and [0.5,1] for a close 
to zero. Define 

+3x 

—3x + 3a + 1 
H a (x) = < +3x — a — 1 
—3a; + 3a + 3 
+3x-2 



< x < \ + ±a, 



| + \a < x < | + |a, 
| + |a<x<g + |a, 

| + \a < x < 1, 



where values are taken modulo 1. Figured shows a graph of H . Let aj e 
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be close to zero, for example (ai, a 2 , 03, 04) = (tc, 2y2, a/3, e)/40. We now construct the map 
Ti from H ai , for i — 1, ... ,4 as follows: 

Ti = H ai (x) 

T 2 = RoH a2 (x) 

T 3 = H^oR- 1 

T 4 = RoH aA oR-\ 

where R : S 1 — )■ 5 1 is i/ie rotation R(x) — x+ 1/4 (mod 1); see FzgweEl 

Let m denote normalised Lebesgue measure on M. To each map T u we associate a 
Perron-Frobenius operator P w : L 1 (M, m) defined by 7^/ = JZy^r^x f(v)/\ DT u (y)\. 
The operator is a linear operator that acts on integrable functions in analogy to the 
action of on points. If / G L l (M,m) represents a density function for an ensemble of 
initial conditions, then V w f represents the density function of the ensemble after the action 
of T w has been applied to the ensemble. The map cocycle $ naturally generates a Perron- 
Frobenius cocycle Vw = V a k-i u o ■ ■ ■ o V a ^ o V u - This composition of k Perron-Frobenius 
operators capture the action on a function / after k iterations of the non-autonomous system. 
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2.2 Continuous time — Flows 



Let F : Q x M — > M. d be a sufficiently regular vector field. More precisely, we suppose that F 
satisfies the conditions of (2j Theorem 2.2.2], which will guarantee the existence of a classical 
solution of the nonautonomous ODE x(t) = F(9 t u,x(t)), t G R. 

To be concrete about the probability space (Q, F, P) in the continuous time setting, we 
may set Q = E C R dl , d\ > 3, and consider an autonomous ODE z = g(z) on S. Denote the 
flow for this ODE by £ : R x H — >■ H and suppose that £ preserves the probability measure 
P; that is, P o £(— t, •) = P for all t G P. Thus, the autonomous, aperiodic flow £ drives the 
nonautonomous ODE 

x(Jt) = F(0 t u,x(t)) = F^{t,z),x(t)). (1) 

We think of points z G H as representing generalised time. We assume that (H, £, P) is 
ergodic in the sense that if £(— i, S) = H for some H C H and for all i > then P(H) = or 
1. 

Denote by : R x S x M ->■ M the flow for ([I]). The flow satisfies f t <fr(t,z,x) = 
F(at,z),<t>(t,z,x)). 

Definition 2. Endow M with the Borel a-algebra and let \x be a probability measure on 
M. We call an invariant measure if /x o <f)(—t, z,-) = \x for all z G E and tel. 

Remark 1. In Definition |2] we are insisting that fi is preserved at all "time instants". As 
in the discrete time setting, more generally one may allow fi = \x z and insist that /j£(_ tjZ ) o 
<f>(—t,z,-) = n z . However, as we will soon begin to focus on coherent sets rather than 
invariant measures, we will restrict the invariant measure to a "time independent" measure 
for clarity of presentation. This is perfectly reasonable for one of the main applications we 
have in mind, namely, aperiodically driven fluid flow where /i = Lebesgue, and volume is 
preserved by the flow at all times. 

Example 2. Consider the following nonautonomous system on a cylinder M = S 1 x [0, 7r]. 
Let ( : R x I 3 4 i 3 denote the flow for the driving system generated by the Lorenz system 
of ODEs with standard parameters a = 10, ft = 8/3, p = 28. 

z\ = (t(z 2 -z 1 )/t (2) 
z 2 = (pzi- z 2 - z 1 z 3 )/t (3) 
4 = {-ftz 3 + ziz 2 )/r. (4) 

It is well known that this Lorenz flow possesses an SBR measure P IJ^j - Let the time- 
dependent vector field F : R x S 1 x [0, it] — > S 1 x [0,n] generate our non- autonomous ODE 
(x(t),y(t))=F(£(t,z),x(t),y(t)). Explicitly, 

x = c — A sin(x — uzi(t)) cos(y) (mod 2n) (5) 
y = A cos(x — uzi(t)) sin(y), (6) 

with c = 0.5, A — 1, v — 0.25. We set initial condition z(0) = (0,1,1.5) and take the 
z\-coordinate of the Lorenz driving system to represent the generalized time for the vector 
field F(£(t, z),x(t),y(t)). We use a scaling factor of r — 6.6685 so that the temporal and 
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Figure 4: Trajectory of the time-dependent system ©-(jS]) driven by the Lorenz system at 
generalized times z) 



spatial variation of z\(t) is similar to that of the "actual" time t. Since F(£(t,z),x,y) is 
differentiable and bounded on M for all t, classical solutions of the nonautonomous ODE 
exist. The system uniquely generates an RDS, see J|J Theorem 2.2.2]. In 

Figure^ we demonstrate a trajectory of three different initial points. 

We may define a family of Perron-Frobenius operators as Vz f(x) = /(</>(— t, z),x)) ■ 



det D(f>(-t, £(t, z),x)\ for t > 0. This family is a semigroup in t as vi tl+t2) f = V^ t ] z) V^'f. 



>(*2) ^(tl) 



3 Galerkin projection and matrix cocycles 

Let B n = sp{xBi '■ Bi E ^3} where 58 = . . . , B n } is a partition of M into connected sets 
of positive Lebesgue measure. Define a projection ir n : L l (M,m) — > B n as 



7T 



Following Ulam 



in the sequel we will consider the finite rank operators rc n V l 



(7) 



(i) 



L \M,m) — > B n and n n Vi 1] : L X (M, m) — » i3 n , and the matrix representations of the re- 
strictions of ir n V^ and to B n . We denote these matrix representations (under mul- 
tiplication on the right) by P(oS) and P(z). Extending Lemma 2.3 in a straightforward 
way to the nonautonomous setting, one has 



and 



P(z) 



ij 



m{Bj n $(-l,au},Bi)) 
~ m(Bi) 

m(£,n(K-U(M,^)) 

m(Bi) 



(9) 



In particular, these matrices are numerically accessible. 
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Remark 2. Note we do not concern ourselves at all with the relationship between Pi and 
VL 1] ; this is a subtle relationship and beyond the scope of this work. See |29j Q2], [TUl E21 
CEH El ?] for work in this direction. 

The matrices P{oS) and P{z) generate matrix cocycles 

:= Pia^u) ■ ■ ■ P(auj) ■ P(u) (10) 

and 

P<*>(z) := Pm - M)) • • • P(£(l, *)) • P(z). (11) 

4 Discretised Oseledets functions and the Multiplica- 
tive Ergodic Theorem 

In periodically driven flows, Liu and Haller [30] and Pikovsky and Popovych [31], observed 
that certain tracer patterns persisted for long times before eventually relaxing to the equilib- 
rium tracer distribution. Pikovsky and Popovych [31] recognised these patterns as graphs of 
eigenfunctions of a Perron-Frobenius operator corresponding to an eigenvalue L < 1. These 
eigenfunctions decay over time and the closer L is to 1, the slower the decay and the more 
slowly an initial tracer distribution will relax to equilibrium. We now develop a framework 
for the considerably more difficult aperiodic setting. 

Consider some suitable Banach space (J 7 , || ■ ||) of real valued functions; J 7 is the function 
class in which we search for slowly decaying functions. Suppose that the norm is chosen so 
that for each oo G f2 and k > 0, the operator V^f 1 is Markov; that is, \\Vu || = 1 for all oo 
and k > 0. For / G J, we calculate the following limit: 

A(w,/)=limBupitog||pW/l|. (12) 

k— »oo ™ 

We refer to X(u, f) < as the Lyapunov exponent of /. If / decays under the action of the 
Perron-Frobenius operators at a geometric rate of r k , < r < 1, then A (a;, /) = logr. The 
closer r is to 1, the slower the decay. The extreme case of r = 1 (no decay) is exhibited when 
/ is the density of the invariant measure \x that is common to all maps in our nonautonomous 
dynamical system. We define the Lyapunov spectrum A(V,u) := {A(cj, /) : / G J 7 }. In the 
aperiodic setting the new mathematical objects that are analogous to strange eigenmodes 
and persistent patterns will be called Oseledets functions. 

Definition 3. Oseledets functions correspond to / for which (i) A (a;, /) is near zero and (ii) 
the value X(u, f) is an isolated point in the Lyapunov spectrum. 

By considering (J 7 , || • ||) = (B n , \\ ■ ||i), the actions of and are described by 
p( k \ui) and P^(z), respectively. We may replace and P { z k) in by P( k \uj) and 
p( k \z), respectively, to obtain a standard setting where the possible values of X(u,f) are 
the Lyapunov exponents of cocycles of n x n matrices, and 

A(P,w) := j Km ilogHP^M/IK = / e B n \ (13) 



S 



and 

A(P,z) := { lim ± log ||P (fc) (z)/||i : / G B n | , (14) 

exist for P almost-all cu G f2, and consist of at most n isolated points, A n < • • ■ < Ai = 0. Of 
particular interest to us is the function /^(w) (or /^(z)) in B n , which represents the function 
that decays at the slowest possible geometric rate A2. 

Remark 3. In certain settings, this matrix cocycle exactly captures all large isolated Lya- 
punov exponents of the operator cocycle V : (BV, || • ||bv) O- One such setting is a map 
cocycle formed by composition of piecewise linear expanding maps with a common Markov 
partition 03 = {B 1} . . . , B n }; see [i~9] . 

The following example illustrates the concept of Lyapunov spectrum and Oseledets func- 
tions in the familiar autonomous setting. For the remainder of this section, we adopt the 
discrete time notation of a and u. 

Example 3 ("Autonomous" single map). In [5] individual maps are constructed for which 
the Perron-Frobenius operator has at least one non-unit isolated eigenvalue when acting on 
the Banach space (BV, \\ ■ \\bv)- A single autonomous map may be regarded as a cocycle 
over a one-point space Q = {u>}, and so we may drop the dependence on u in notation. 
Keller [27^ shows that for a piecewise expanding map T of the interval I , the spectrum of 
the associated Perron-Frobenius operator V has an essential spectral radius p eS s(P) equal to 

I h I 1/^ 

the asymptotic local expansion rate sup xg/ lim^oo 1/DT (x)\ , and that there are at most 
countably many spectral points, each isolated, of modulus greater than p ess (V). In order to 
have an isolated spectral point, we construct a map of S 1 which has an almost-invariant set 
(see Definition^). The relation between almost-invariant sets and isolated eigenvalues was 
noted in Q/. Consider the partition 03 = {B t : i — 1, . . . ,6}, where Bi = ({i — l)/6,i/6). 
Given a G Z 6 , any map T : S* 1 — >■ S 1 defined by 

T(x) = 3x- (i- l)/2 + a</6 (mod 1), x G B { (15) 

is Markov with respect to 03. Here we take a = (0,0,1,4,3,3); see Figure\5[ Notice that 
there is a low transfer of mass between the two intervals [0, 1/2] and [1/2, 1]. Since 23 is a 
Markov partition for T, the space of characteristic functions Bq = {xBi '■ i — 1, • • • > 6} is 
an invariant subspace o/BV for the Perron-Frobenius operator V ofT. Thus the action of 
Vuj = V on T = Bq is represented by the matrix 



P = P(u) 



( 1 


1 





1 \ 


1 


1 


1 





1 


1 


1 











1 


1 1 











1 1 1 


\o 








1 1 l) 


hL 2 




(1- 


V V2)/3,L 



(16) 



which has non-zero eigenvalues L\ = 1,— , . v _,, , ■ v .. . .. , 

is piecewise affine with constant slope 3 and so the logarithm of the local expansion rate is 
log(l/3). 
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Figure 5: Graph of T and Oseledets function f 2 . 



The eigenvalue L 2 ~ 0.805 of P thus gives rise to an isolated point X 2 ~ log 0.805 in 
the Lyapunov spectrum A(V). The corresponding Oseledets function f 2 is given by f 2 {x) = 
Y^h=i w 2,iXBi( x ), where w 2 is the eigenvector of P corresponding to the eigenvalue L 2 , see 
Figured Since \L$\ ~ 0.138 < 1/3, this means that logL 2 is the unique isolated Lyapunov 
exponent in A{V). Note that the set {f 2 > 0} corresponds to the set [0, 1/2]. We will discuss 
this property further in Section [7| 



Example 4 (Periodic map cocycle). We construct a periodic map cocycle from a collection 
of maps with a common Markov partition. The map cocycle is formed by cyclically composing 
three maps of S 1 . Consider the sequence space Q = {u £ |l,2,3| z : Vi G Z, 
(0 1 0\ 

where M = I 1 . We consider T = {Tj : j 
\1 Qj 

with parameter where 



{u e {1,2,3} Z : G 1,M m+1 = 1} 
1,2,3}, where Tj is given by < f73]) 



a« = (3,2,2,0,5,5), 



,( 2 ) 



(2,1,4,5,4,1) 



,(3) 



;i, 3, 3, 4, 0,0), 



see Figure® As in Example^ we look for Lyapunov exponents that are strictly greater than 
the logarithm of the asymptotic local expansion rate 



sup lim \1/D(T 3 o T 2 o T\ 



,fc|l/3fc 



(17) 



As each map Tj is piecewise afftne with constant slope 3, the logarithm of the local expansion 
rate is log(l/3). Note also that T\ approximately maps [0,1/2] to [1/3,5/6], T 2 then maps 
[1/3,5/6] approximately to [0,1/6] U [2/3,1], and finally T3 maps [2/3,1/3] approximately 
back to [0, 1/2]. Each map Tj leaves the space B e from Example invariant, and thus the 
Perron- Frobenius operator Vj ofTj restricted to Bq has matrix representation Pj, where 3Pj, 
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Figure 6: Graphs of T\, T 2 and T 3 . 



j = 1,2,3, are respectively 
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1 







1 





1 
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1 
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1 




1 


1 


1 











1 


1 


1 













1 





1 





1 










1 


1 


1 








V 1 











1 


1 1 









1 


1 


1 


0/ 




\o 


1 


1 


1 





0/ 



The triple product p( 3 \tu) = P 3 P 2 P\ has non-zero eigenvalues Li — 1, L 2 — (13 + v / 233)/54 
and L 3 = (13— v^233)/54. Since L 2 ~ 0.523, its associated eigenvector w 2 satisfies X(u , w 2 ) = 
log v^2 > log(l/3). Since P^\a k oj), k = 1,2, are cyclic permutations of the factors of 
p( 3 \cu), they share the same eigenvalues, and in particular L 2 . Thus (1/3) log L 2 is an iso- 
lated Lyapunov exponent of A(V,u) for each u G Q. Associated to the eigenvalue L 2 , the 
matrices P(a k u), k = 0,1,2, have corresponding eigenvectors w 2 (a k u). The three vectors 
w 2 {a k u), k = 0, 1, 2, generate the periodic Oseledets functions f 2 (cr k uj) = Y^h=i u '2,i(o" fc ^ mod 
see Figure\7\ Note that the sets {f 2 (a k uj) > 0}fc=o,i,2 correspond to the sets [0, 1/2], [1/3, 5/6], 



f 2 {au) 



f 2 (a 2 u) 



Figure 7: Oseledets functions f 2 (o~ k u) for k — 0, 1, 2. 

and [0, 1/6] U [2/3, 1], respectively. We will discuss this property further in Section^ Tor 
another such example see \Wj - See also WD^ for a detailed example of similar calculations 
for a periodically driven flow. 
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In the nonautonomous setting, we can no longer easily construct Oseledets functions as 
eigenf unctions of a single operator, or eigenvectors of a single matrix. In fact, the Oseledets 
functions are themselves (aperiodically) time dependent in the nonautonomous setting. Our 
model of Oseledets functions for nonautonomous systems is, as the name suggests, built 
around the Multiplicative Ergodic Theorem, see e.g. [2, Chapter 3, §4]. We now state a 
strengthened version [19] of the Multiplicative Ergodic Theorem that we require for our 
current purposes. 

Theorem 1 ([19]). Let a be an invertible ergodic measure-preserving transformation of the 
space (Q,7i,¥). Let P : Q — > M n (R) be a measurable family of matrices satisfying 

J log + ||P(w)|| dP(w) < oo. 

Then there exist Ai > A 2 > • • ■ > A^ > — oo and dimensions mi, ... , m t , with m x + ■ • - + m^ = 
n, and a measurable family of subspaces Wj(u) C R™ such that for almost every u G Q the 
following hold: 

1. dimWj(co) = rrij; 

2. R n = e- =1 W»; 

3. P(u)Wj(uj) C Wj(au) (with equality if \j > -oo); 
4- for all v G Wj(co) \ {0}, one has 

lim (1/jfe) log WPia^uj) ■ ■■P(au) ■ P(u)v\\ = Xj. 

k—>oo 

The subspaces Wj(u>) are the general time-dependent analogues of the vectors w 2 and 
W2{cr k u), k = 0, 1, 2 of Examples [3] and HI respectively. We may explicitly construct a slowest 
decaying discrete Oseledets function as fz{w) '■= Y^Ji=i w 2,i(^)XBi, where itf 2 (w) G W^u;). In 
the sequel, for brevity we will often call Wj (u) a subspace or a function, recognising its dual 
roles. 

Remark 4. We remark that if £ > 2, m 2 = 1, and A 2 > — oo, the family of vectors 
{/ 2 (cr fc u;)}fc>o is the nnzgw^l (up to scalar multiples) family of vectors in B n with the prop- 
erties that 

1. hm fc ,_ >00 (l/A;')log||P( fc ')(a;)/ 2 (a fc a;)|| 1 = A 2 , k > 0, 

2. ^nV LU f 2 {o- k uj) = a k f 2 (o- k+1 uj) for some a k ^ 0, k > 0. 

2 Assume there is another family {w 2 (cr k uj)}k>o ^ {w2(c k ou)}k>o (up to scalar multiples) with these 
properties. Then w' 2 (cr k uj) — Ylj=2 a k,jWi(a k uj) for some a k ,j, j = 2, . . . ,£, with a k _ 2 ^ 0. WLOG assume 
®k,2,®k,j' 7^ for some 2 < j' < £ and all k > 0, but that a k ,j = for all j ^ 2,j' and all k > 0. Then 
77i 2 > 2 in Theorem [TJ a contradiction. 
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Remark 5. Theorem [I] strengthens the standard version of the MET for one-sided time with 
noninvertible matrices (see e.g. [21 Theorem 3.4.1]) to obtain the conclusions of the two-sided 
time MET with invertible matrices (see e.g. [2j Theorem 3.4.11]). In [21 Theorem 3.4.1], the 
existence of only a flag W" 1 = Vi(co) D • ■ ■ D Vi{u) of Oseledets subspaces is guaranteed, while 
in [21 Theorem 3.4.11], the existence of a splitting W\{oS) © ■ • • © Wi(u) = W 1 is guaranteed. 
Theorem [1] above demonstrates existence of an Oseledets splitting for two-sided time with 
noninvertible matrices. This is particularly important for our intended application as the 
projected Perron-Frobenius operator matrices are non-invertible. Recent further extensions 
[?] prove existence and uniqueness of Oseledets subpsaces for cocycles of Lasota-Yorke maps. 



5 Numerical approximation of Oseledets functions 

In the autonomous and periodic settings we have seen in Examples [3]andH]that the subspaces 
W 2 (^) = sp{w 2 (^)} were one-dimensional, and that the vectors 102(00) could be simply 
determined as eigenvectors of matrices. For truly nonautonomous systems (those that are 
aperiodically driven), the Oseledets splittings are difficult to compute. In this section we 
outline a numerical algorithm to approximate the Wj(uo) subspaces from Theorem [TJ The 
algorithm is based on the push- forward limit argument developed in the proof of Theorem [H 
To streamline notation, we describe the discrete time and continuous time setting separately. 

5.1 Discrete time 

We first describe a simple and efficient method to construct the matrix P(uo) defined in (JHJ). 
Algorithm 1 (Approximation of P(o~~ k uo)ij, < k < N). 

1. Partition the state space M into a collection of connected sets {-Bi, . . . , B n } of small 
diameter. 

2. Fix i, j , and k and create a set of Q test points Xj tl , . . . , x^q G Bj that are uniformly 
distributed over Bj . 

3. For each q = 1, . . . , Q calculate yj iq = T c -h^Xj A . 

4. Set 

pr„-KA #{g : y j>q e B,} 

P(a uo)ij = (19) 

We now describe how to use the matrices P(to) to approximate the subspaces Wj(co). 
An intuitive description of the ideas behind Algorithm [2] immediately follows the algorithm 
statement. 

Algorithm 2 (Approximation of Oseledets subspaces Wj(oo) at to G O.). 

1. Construct the Ulam matrices P( m \o~ n uj) and P( N \a~ N 00) from / fTPj) and ^^) for 
suitable M and N. The number M represents the number of iterates over which one 
measures the decay, while the number N represents how many iterates the resulting 
"initial vectors" are pushed forward to better approximate elements of the Wj(uo). 
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2. Form 

¥ M \a- N uj) := (P^ M \a- N uj) T P^ M \a- N uj)) l ' 2M 
as an approximation to the standard limiting matrix 

B(a- N u) : = lim (P^ (a' N u) T P^ M \a~ N u)) 1/2M 

appearing in the Multiplicative Ergodic Theorem (see e.g. Jlj, Theorem 3.4-l(i)])- 

3. Calculate the orthonormal eigenspace decomposition of^>( M \a~ N uj), denoted by Uj(a' 
j = !,...,£. We are particularly interested in low values of j , corresponding to large 
eigenvalues Lj. 

4- Define Wj (u) := P^ N ^ (o~~ N uj)Uj M \a~ N u) via the push forward under the matrix 
cocycle. 

(u)) is our numerical approximation to Wj(u). 

Here is the idea behind the above algorithm. If we choose M large enough, the eigenspace 
Uj M \a~ N ui) should be close to the limiting (M — > oo) eigenspace Uj(a~ N uj). Vectors in the 
eigenspace Uj(a~ N u) experience stretching at a rate close to Lj. Note that the eigenspace 
Uj M \a~ N u) is the j th singular vector of the matrix P^ MS> (a^ N u), which experiences a "per 
unit time" average stretching from time — N to —N + M of Lj. Choose some arbitrary v G 
Uj{a~ N uj) and write v = Y^j'=j w j' with uy G Wf(a~ N u). Pushing forward by p( N \a~ N uj) 
for large enough iV will result in \\P( N \a~ N uj)wj\\ dominating \\P( N \a~ N uj)wji\\ for j < j' < 
£. Thus, for large M and iV we expect Wj M ' N \u) to be close to Wj{uS). 

Remarks 1. 

1. Theorem □ states that W^' N \u) ->■ Wj(u) as N ->■ oo. 

2. This method may also be used to calculate the Oseledets subspaces for two-sided linear 
cocycles, and may be more convenient, especially for large n, than the standard method 
of intersecting the relevant subspaces of flags of the forward and backward cocycles. 

The numerical approximation of the Oseledets subspaces has been considered by a variety 
of authors in the context of (usually invertible) nonlinear differentiable dynamical systems, 
where the linear cocycle is generated by Jacobian matrices concatenated along trajectories of 
the nonlinear system. Froyland et al. [18] approximate the Oseledets subspaces in invertible 
two-dimensional systems by multiplying a randomly chosen vector by (o~~ n uj) (pushing 
forward) or P^~ N \a N u) (pulling back, where P^ N \a N u) = p-\u) ■ ■ ■ P- 1 (a N ~ 1 u). Tre- 
visan and Pancotti [?] calculate eigenvectors of ^ M ^{u) for the three-dimensional Lorenz 
flow, increasing M until numerical convergence of the eigenvectors is observed. Ershov and 
Potapov [12] use an approach similar to ours, combining eigenvectors of a ^f^ M ' with pushing 
forward under P^ N \ Ginelli et al. [23] embed the approach of [18] in a QR-decomposition 
methodology to estimate the Oseledets vectors in higher dimensions. In the numerical ex- 
periments that follow, we have found our approach to work very well, with fast convergence 
in terms of both M and N. 
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5.2 Continuous time 



As our practical computations are necessarily over finite time intervals, from now on, when 
dealing with continuous time systems, we will compute P^(z) as it n Vz rather than as 
P(£(fc — 1, z)) ■ ■ ■ P(£(l, z)) ■ P{z). If the computation of n n V^ can be done accurately (this 
will be discussed further in Section |53|) . then this representation should be closer to Vz as 
there are fewer applications of 7r n . 

We first describe a simple and efficient method to construct the matrix P(co) defined in 

®. 

Algorithm 3 (Approximation of p( M \£(-N, z)), N > 0). 

1. Partition the state space M into a collection of connected sets {B\, . . . , B n } of small 
diameter. 

2. Fix i,j , and z and create a set of Q test points Xj t i, . . . , x^q G Bj that are uniformly 
distributed over Bj . 

3. For each q = 1, . . . , Q calculate yj A = (f)(M,£(—N, z),Xj iq ). 

4. Set 

pW(e(-A,z)) tJ = * {q ■ Vi « E Bi} (20) 

The flow time M should be chosen long enough so that most test points leave their 
partition set of origin, otherwise at the resolution given by the partition {£>!, . . . , B n }, the 
matrix p( M \£(— N, z)) matrix will be too close to the nx n identity matrix. If the action 
of (j) separates nearby points, as is the case for chaotic systems, clearly the longer the flow 
duration M, the greater Q should be in order to maintain a good representation of the 
images <fi(M,£(—N, z), Bj) by the test points. 

Algorithm 4 (Approximation of Oseledets subspaces Wj(z) at z G 5.). 

1. Construct the Ulam matrices P^ M \^(—N, z)) and P^ N \^(—N, z)) from (E?J) for suit- 
able M and N . The number M represents the flow duration over which rate of decay is 
measured, while the number N represents the duration over which the resulting "initial 
vectors" are pushed forward to better approximate elements of the Wj(z). 

2. Form 

¥ M \a-N,z)) : = (p( M \a-N,z)) T p^\a-N,z))) 1/2M 



as an approximation to the standard limiting matrix 

B(£(-N,z)) := lim (P^\^-N, z)) T P^\^~N, z))f 2M 
appearing in the Multiplicative Ergodic Theorem (see e.g. |H Theorem 3.4- l(i)])- 



3. Calculate the orthonormal eigenspace decomposition of \I/( M )(£(— A, z)) , denoted by 
Uj M \^(— A, z)), j — 1, . . .,£. We are particularly interested in low values of j, corre- 
sponding to large eigenvalues Lj. 
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4- Define Wj M (z) := P( N \£,(—N, z))Uj (£(—N, z)) via the push forward under the 
matrix cocycle. 

5. W 3 iM ' N) (z) is our numerical approximation to Wj(z). 

5.3 Continuity of the Oseledets subspaces in continuous time 

When treating continuous time systems, one may ask about the continuity properties of 
Wj' N \z) in z. In the following we suppose that W2, M ' N \z) is one-dimensional. For large 

M and N, W^' 1 ^ (z) will approximate the most dominant Oseledets subspace at time z. 
Suppose that we are interested in how this subspace changes from time z to time £ (8, z) 
for small 8 > 0. There are two ways to obtain information at time £(8,z). Firstly, we can 
simply push forward W^'^iz) slightly longer to produce W% ) M,N z)). Secondly, we 

can compute \I/( M ) slightly later at time £(8,z) to produce W2 M ' N \£(8, z)). 

To compare the closeness of W2 (z) to z)) and W% i z)), we 

represent each as a function and make a comparison in the L 1 norm. We assume that 
u \ \(£(- N , z )) is one-dimensional and define f n ,t(-N,z),M = Th=i( u 2^(- N ^ z )))iXB, G 
L\M,m) where uf\^{-N,z)) G U^ M) ^{-N, z)) is scaled so that \\f n ,t(-N,z),Mh = 1 - 
Let f n ,z,M,N = n n vj£} N z) fn,t(-N,z),M- Note that fn,z,M,N = Yh=i (w^ ( z ))iXB, for some 
wi M ' N \z)eW^ M ' N) (z).[ 

We firstly compare f n ^(s,z),M,N+s and f n>z>M ,N- 

Proposition 1. \\f n ,as,z),M,N+5 ~ /n,z,M,Jv||i ->■ as 8 -)> 0. 

Proof. Note that fn£{5,z),M,N+5 = ^nP^^fn^-N ~,z),M While f n ,z,N,M = ^nP^} NtZ )fn^(-N,z),M- 

The proof will follow from the result that P^ is a continuous semigroup; that is, lim^o \\Vl f- 
/||i = for all t G M, / G L 1 (M, m). 

Lemma 1. \\vi 5) f - f\\i -> as 8 ^ for all z e E and f e L 1 . 

Proof. The proof runs as a non-autonomous version of the discussion in Remark 7.6.2 [28J. 

Note that T>jPf(x) = f(<f>(-5,£(5,z),x)) ■ detD<f>(-5,£(5,z),x), where <f>(-5,£(5,z),-) de- 
notes the flow from £ (5, z) in reverse time for duration 8. For the moment consider continuous 
/. Since x \- > <p(s, z, x) is at least C 1 for each s, z (the derivative of wrt to x is continuous 
with respect to s and x for each fixed z) by [21 Theorem 2.2.2 (iv)] and M is compact, 
Pz^ f(x) — > /(x) uniformly in x as 8 — > 0. Thus \\Vz^ f — f\\i — > as 8 — )■ 0. Since the 
continuous functions are dense in L p , 1 < p < oo as M is compact (see e.g. [H] Lemma 
IV. 8. 19), one can L 1 approximate any L l f by a continuous function and thus the result 
holds for all L 1 functions /. □ 

Thus the result follows using Lemma [1] and the fact that ||7r n ||i = 1. □ 

Now, let's compare f n ,^s,z),M,N and f n , z ,M,N- 

Proposition 2. \\f n ,t(s,z),M,N ~ fn,z,M,N\U as 8 ^ 0. 
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Proof. This result is more difficult to demonstrate as we need to firstly compare \I/( M ) (£(— N, z)) 
with ty( M > (£(—N + S, z)). To this end, consider 

< llir.ll, (ll^fW " + IKS,/ " <L,/H 



The right hand side converges to zero as 5 — > by Lemma [TJ This result implies that 
||p(M)^_jy^ ^ _ p( M )(£(_j\r _|_ ^ z )) || —7- as 5 —7- in whatever matrix norm we choose. 
Thus \\^ M \^-N,z)) 2M -^ M \^(-N + S,z)) 2M \\ = \\P( M \Z(-N,Z)) T (P( M \Z(-N,z)) - 
p(M)(£(-N + 5, z))) + (P( A/ )(^(-iV, z)) T - p( M )(e(-iV + 5, z)) T )p( M )(£(-iV + 5, z))\\ ->■ 
as 5 —7- 0. By standard perturbation results, see e.g. [26, Theorem II.5.1], this implies that 
eigenvectors U2 (z) and (£(<5, z )) are close for sufficiently small 5. Thus f n ,£(-N,z),M 
and fn,t(-N+s,z),M are close in L 1 norm. Now we need to push both of these forward by 
7i n V(^(—N, z)Y N \ This will not increase the norm of the difference at all, so \\f n ,£(5,z),M,N ~ 
fn,z,M,N\\i will also be small. □ 



5.4 Oseledets functions for a ID discrete time nonautonomous 
system 

We now examine the Oseledets functions for the system defined in Example [TJ We consider 
the approximation 7TiooPo; of rank 100, which we obtain by Galerkin projection. We denote 
by P(oo) G M 100 x M 100 the Ulam matrix representing the action of 7TiooPw on functions 
/ G B wo := sp{x[(i-i)/ioo,i/ioo)j* = !,•••, 100}. The matrices P((j- k uj),k = -10, . . . , 10 are 
constructed by following Algorithm 1 using Q = 100. 

We look for Oseledets functions for a particular aperiodic sequence u. To generate an 
aperiodic sequence, let r G {0, 1} N be the binary expansion of l/y/3. Extend r to an element 
of {0, 1} Z by setting Tj = for all i < 0. Define Co>j_ 2 5 = 1 + 2rj + r i+1 for each i. Then oj G Q 
and the central 21 terms of u are 

u = (...,2,3,1,2,4,4,3,2,3,1,1,2,3,2,4,3,1,2,3,1,1...), (21) 

where the dot denotes the zeroth term coq = 1. 

We calculate the eigenvalues of (F^irtj^Pl 20 '^ 1 ^)) 1 / 40 , where p( 20 \o- w cu) is 
defined as in (JTUJ), and find the top three to be 

L x « 1.00, La « 0.84, L 3 « 0.46. 

As in Examples [3] and HJ the maps Tj are piecewise affine with constant slope three, and so 
p(uj) = 1/3. Thus logL 2 and logL 3 may approximate isolated Lyapunov exponents in A(V). 

We follow Algorithm 2 to approximate the second Oseledets subspace W% (a k u)) for 
k — 0, ... ,5, using (M, N) = (20, 10), see Figure [HJ In order to confirm the effectiveness of 
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Figure 8: The Oseledets function approximations f^ 1 {a k uS) for M = 20, N = 10, and 
k — 0, ... ,5, along with optimal thresholds (shown in dashed green), see Section 17.21 

A(iV) 




N 



Figure 9: A graph showing A(N) for N = 1, ... ,19. 
Algorithm 2 we calculate the L 1 distance A(iV) between the normalisations of the vectors 



w 



(2N,N) , 



au) and P(u>)w^ N ' N \u>), for N = 2, . . . , 19 with M = 40. By property 3 of Theorem 
[T]this distance should be small if the family W^uS) is well approximated. A logarithmic plot 
of A (AT) against N, see Figured shows the fast convergence of wf N,N \u) to an Oseledets 
subspace. In Section 17.21 we will see how to extract coherent sets from these functions. 



5.5 Oseledets functions in a 2D continuous time nonautonomous 
system 

We consider the following nonautonomous system on M = [0, 2ir] x [0,7r], t G M + : 



x = c — A sin(x — vt) cos(y) (mod 2ir) 
y = A cos(x — vt) sin(?/) 



(22) 
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This equation describes a travelling wave in a stationary frame of reference with rigid bound- 
aries at y — and y = tt, where the normal flow vanishes [33] 139] . The streamfunction 
(Hamiltonian) of this system is given by 

s(x, y, t) = —cy + A sin(x — vt) sin(?/). (23) 

We set c = 0.5, A = 1, and the phase speed to v = 0.25. The velocity field is 27r-periodic 
in the x-direction, which allow us to study the flow on a cylinder. The velocity fields in 
a comoving frame for these parameters are shown in Figure [TDJ The closed recirculation 
regions adjacent to the walls (y = and y = tt) move in the positive x-direction and are 
separated from the jet flowing regime by the heteroclinic loops of fixed points, which are 
given below. 

This model can be simplified to an autonomous system with a steady streamfunction 
in the comoving frame by setting X = x — vt and Y = y. The steady streamfunction 
is then given by S(X,Y,t) = -(c - v)Y + Asm(X) sin(F). Let X s = sin~ 1 ((c - v)/A) 
and Y s = cos _1 ((c — v)/A). In the comoving frame, the recirculation region at the wall 
Y = contains an elliptic point q± = (ir/2,Y s ) and is bounded by the heteroclinic loop 
of the hyperbolic fixed points p\ = (X s , 0) and p 2 = (tt — X s ,0). Similarly, those elliptic 
and hyperbolic points at the wall Y = tc are g 2 = (37r/2,7r — Y s ), p 3 = (n + X s ,n), and 
P4 = (27r — X s ,7r), respectively, see Figure [TOl One may observe that there is a continuous 
family of invariant sets in the comoving frame as any fixed level set of the streamfunction 
bounds an invariant set. In a stationary frame these elliptic and hyperbolic points (and 
their heteroclinic loops) are just translated in the x— direction. That is, any fixed level set of 
the time- dependent streamfunction f)23p is a (time-dependent) invariant manifold. We note, 
however, that the recirculation regions are distinguished from the remainder of the cylinder 
as they are separated from the jet flowing region, which has a different dynamical fate. In 
the subsequent sections we will perturb this somewhat "degenerate" system to destroy the 
continuum of invariant sets in the comoving frame and produce a small number of almost- 
invariant sets (see Definition [5]) in the comoving frame, or coherent sets in the stationary 
frame) . 



5.5.1 A coherent family: Mixing case 

We modify the traveling wave model in the previous section to allow mixing in the jet flowing 
region. We add a perturbation to the system in the following way: 

x = c - A(z(t)) sin(x - uz{t)) cos(y) + eG(g(x, y, z{t))) sin(5(t)/2) 
y = A(z(t))) cos(x — vz{t))) sm{y). 

Here, z(t) = 6. 6685^1 (t), where z\(t) is generated by the Lorenz flow in Example [2] with 
initial point z(0) = (0,1,1.5), A(z(t)) = 1 + 0.125 sin(V5£(t)), G(ip) := l/(^ 2 + l) 2 and 
the parameter function ip = g(x,y,z(t)) := sin(a; — uz(t)) sin(y) + y/2 — ir/4 vanishes at 
the level set of the streamfunction of the unperturbed flow at instantaneous time t=0, i.e., 
s(x,y,0) = 7r/4, which divides the phase space in half. We set e — 1 as this value is 
sufficiently large to ensure no KAM tori remain in the jet regime, but sufficiently small to 
maintain islands originating from the nested periodic orbits around the elliptic points of the 
unperturbed system. 
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Figure 10: Vector fields in the comoving frame for the travelling wave flow (|23l . for A = 
1.0 and c = 0.5. The red dots are the hyperbolic fixed points that are connected by the 
heteroclinic loops. The blue dots are elliptic points in the centre of recirculation regions. 

We applied Algorithm H with n = 28800, M = 80, N = 40, z = (0, 1, 1.5), and Algorithm 
H]for z = (0, 1, 1.5) and z = £(10, (0, 1, 1.5)). By using a relatively large number of test points 
per grid box (n = 400 points per box Bj) we are able to flow for M = 80 units of time and still 
well represent 0(80, £(—40, z), Bj). Figure ITT1 shows that the resulting Oseledets functions 
highlight the remaining islands in the perturbed time-dependent flow. We calculate the 
eigenvalues of (P( 80 )(£(-40, z)) T o p( 80 )(£(-40, z)) 1 / 2 , where P( 80 )(£(-40, z)) is defined as in 
( |20~1) . and find the top three to be 

Li « 1.1100, L 2 « 0.9691, L 3 w 0.9676. 

By part 3 of Theorem [1] (bundle invariance of W 2 (z)) we should have P( 1Q \z)W% f° (z) ~ 
W2 8 °' i0 \£,(W, z)). This is demonstrated in Figure [TT1 by comparing subplots (e) and (f). In 
Section 17.31 we will see how to extract coherent sets from these Oseledets functions. 

6 Invariant Sets, Almost-Invariant Sets, and Coherent 
Sets 

We begin by briefly recounting some of the background relevant to almost-invariant sets. If 
$ (resp. 0) is autonomous, then Q (resp. S) consists of a single point, and we may write 
$(— l,co,x) = $(— l,x) (resp. <p(—t,z,x) = <fi(—t,x)). 

Definition 4. In the autonomous setting, we call A an invariant set if $(—1, A) = A (resp. 
<j)(-t,A) = A for all t > 0). 

The following definition generalises invariant sets to almost-invariant sets. In the con- 
tinuous time case we define: 

Definition 5. Let [i be preserved by the autonomous flow <fi. We will say that a set A C M 
is po-almost-invariant over the interval [0,r] if 
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Figure 11: (a) Graph of approximate Oseledets function w!f 0,m \z) produced by Algorithm 



HI (b)-(e) Pushforwards of W< 



(80,40) . 



W% ; 40 ^(£(10, z)) produced independently by Algorithm HI compare with (e) 



via multiplication by P^ T \z) for r = 2.5, 5, 7.5, 10. (f) 



P»M) := ^ > po 



(25) 



for all s G [0,r], 
2. A is connected. 

If A C M is almost-invariant over the interval [0, r], then for each s G [0, r], the proba- 
bility (according to /i) of a trajectory leaving A at some time in [0, s], and not returning to 
A at time s is relatively small. In the discrete time setting, r = 1, and the obvious changes 
are made in Definition [5j By convention we ask that A is connected; if A is not connected, 
we consider each connected component to be an almost-invariant set for suitable po- 

We now begin to discuss the nonautonomous setting. The notion of an invariant set is 
extended to an invariant family. 

Definition 6. 

1. Discrete time: We will call a family of sets {A^^}, A a k u C M, u G Q, k G Z an 

invariant family if <f>(—k, u, AJ) = A a -k w for all u G £1 and k G Z + . 

2. Continuous time: We will call a family of sets {^(t,«)}, -^(t,*) C M, z G 5, £ G R 



an invariant family if «/>(—£, -2, = A. 



for all z G H and t G 



21 



Motivated by a model of fluid flow, we imagine coherent sets as a family of connected 
sets with the property that the set is approximately mapped onto A^k^ by k iterations 
of the cocycle from "time" oj\ that is, u, A w ) ~ A a h w . The definition of coherent sets 
combines the properties of almost-invariant sets and an invariant family. As we now have 
a family of sets we require one more property beyond those of Definition \5\ in addition to 
modifying the almost-invariance property. In the continuous time case we define: 

Definition 7. Let p be preserved by a flow and < po < 1. Fix a z G S. We will say 
that a family {^(t z )}t>o, A^ iz ) C M, t > is a family of po- coherent sets over the interval 
[0,r]if: 

1. 

P/A A £(t,z)> •= 7-; v > Po, (2oJ 



for all s e [0, r] and t > 0, 

2. Each t > is connected, 

3. /i(v4 5(M) ) = p(A^ V)Z )) for all t,f > 0, 
In discrete time, we replace ( 1261) with 



P„(A0 := ^ > Po, (27) 

r necessarily becomes 1, and we make the obvious changes to the other items in Definition 

m 

We remark that by selecting some A C M of positive p measure and defining A^n z ) : = 
(pit, z, A), t > 0, the family {A^ t ^} t > , is a family of 1-coherent sets. Such a family is not of 
much dynamical interest, as there is nothing distinguishing this family from one constructed 
with another connected subset A' C M. We are not interested in these constructions of 
coherent sets, and in practice the numerical algorithm we present in the next section is 
unlikely to find such sets for chaotic systems. 

7 Coherent sets from Oseledets functions 

We wish to find a family of sets {A z } so that 

(a a \ ._ M4n^(-s,e(s,z),%, 2) )) 

Pp(A x ,A£ M ) .- j^—t (28) 

is large for s e [0, r]. We may rewrite the RHS of ( 1281) as 

Xa z ■ Xct>(-BA e( s,,)) d P ) IvW = ( / V'f ) x Az ■ Xa^ z) dp ) /p(A z ). (29) 



For (129]) to be large we require VzXa z ~ Xa 
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Let us now make a connection with the Oseledets functions f 2 (z) = Yl7=i w 'iA z )XB i 
where w 2 (z) G W 2 (z) obtained in Algorithm HJ In the following discussion, we scale 
f 2 (z) so that ||/2(^)||i = 1 for all z G S. To convert the family of Oseledets functions into 
a family of coherent sets, we modify a heuristic due to [7J that has been successfully used 
in the autonomous setting. The heuristic is to set A z = {f 2 (z) > 0}, z G H. We show that 
Vz f 2 (z) — / 2 + (£(s, z )) i s small for moderate s and large A 2 ; we then heuristically infer that 
V { z s) Xa z ~X% s ,* r 

Proposition 3. Let X 2 = lim s ^. 00 (l / s) log ||7- , i s ' > /2 (^) || < be the second largest Lyapunov 
exponent from Theorem^ and f 2 (z) G W 2 (z) a corresponding Oseledets function, normalised 
so that ||/2(z) ||i = 1. Given an e > there is an S > so that s > S implies \{Pz f 2 { z ) — 



Proof. Given e > we know that there exists S > such that for all s > S one has e 2 e < 



\\rrf2{z)t' s < 1. Since V { z s) h{z) = {V { z s) f 2 (z)) + - {Vf h{z)Y and / VT f 2 (z) dm = 
0, one has ||Pl s) / 2 WI|i = J\V ( z s) } 2 {z)) + + {V^ h{z))~ dm = 2 j (V [ z s) f 2 (z)) + dm. Thus 
J(V ( z s) f 2 (z)) + dm > e^^ s /2. Since (V [ z s) f 2 (z)) + < V ( z s) f 2 + (z), one has \\V ( z s) f 2 + (z) - 
(Vi s) f 2 (z))+W = jV ( z s) f 2 + {z) - (V { z s) f 2 (z)) + dm. As \\f 2 (z)\\ = 1 and j f 2 (z) dm = 0, we 



have / f 2 (z)dm = 1/2 and since V z preserves integrals, J V z f 2 (z)dm = 1/2. Thus, 



The preceding discussion heuristically addresses item 1. of Definition [71 Regarding item 
2 of Definition [7J as we are extracting the sets A z from the Oseledets functions f 2 (z), the 
connectivity of the sets will depend on the regularity of the Oseledets functions. This is a 
delicate question and relatively little can be said formally at present. In the autonomous 
case, roughly speaking, one expects smooth eigenfunctions for Perron-Frobenius operators of 
smooth expanding systems [27_IES], and eigendistributions (smooth in expanding directions, 
distributions in contracting directions) in uniformly hyperbolic settings [3]. These properties 
may carry over to the non-autonomous setting; recent results in the bounded variation 
setting show they do [?] . If a small amount of noise is added by postmultiplying the Perron- 
Frobenius operator by a smoothing (e.g. diffusion) operator, then the Oseledets functions 
must be smooth. This physical addition of a small amount of noise is one way to guarantee 
regularity of the Oseledets functions and connectivity of the associated coherent sets. 

Finally we note that if \x = m one has J f 2 (z)(x) d[x{x) = and so we must have n(A z ) = 
1/2 for all z G 5. Thus, item 3. of Definition [7J is satisfied by the choice A z = {f 2 (z) > 0}. 
If fi 7^ m, then it may be necessary to further tweak the choice of the A z to ensure that item 
3. of Definition [7J is satisfied. This additional tweak is described in Algorithm 

7.1 A numerical algorithm 

For a fixed time z G 5, we seek to approximate a pair of sets A z and A^ TtZ ) for which 



/ 2 + (£M))Hi<(l-e( A2 -^)/2. 




jn s) f 2 (z) - {Vrftiz)) dm < (1 - e^>)/2. 



□ 



Pn{A z , A^ T>z) ) : 



/j(4n^(-r,^(r,2),%, 2) )) 



(30) 
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is maximal. The quantity p^(y4 2 , A^ T)Z )) is simply the fraction of p-measure of A z that is 
covered by a pullback of the set A^ T ^ over a duration of r. For maximal coherence, we 
wish to find pairs A z , A^r T>z \ that maximise p^(A z , A^ T>Z )). We present a heuristic to find 
such a pair of sets based upon the vectors W^ M ' N \z) and W^ M,N ^ (£(r, 2)) corresponding to 
some Lyapunov spectral value A close to 0. This heuristic is a modification of heuristics to 
determine maximal almost-invariant sets, see [171 [T5| [20] . In the terminology of the prior 
discussion in §7, rather than setting A z := {f(z) > 0}, we allow A z := {f(z) > c} or 
A z := {f(z) < c} for some c G K in the hope of finding A z , A^ TiZ \ with an even greater value 
of Pn(A z , A^/ T>z \). This additional flexibility also permits a matching of fi(A z ) and fi{A^f TjZ \). 

Algorithm 5 (To determine a pair of maximally coherent sets at times z, £(t, z)). 

1. Determine W^ M ' N \z) and W^ m,n \^(t, z)) for some r > according to Algorithm^ 

2. Set A+(c) = Ui:w/(Af,iv) (2)>c 5i and A+ (r>2) (c) = Ui:w(«. N )(£(r, z ))>c 5 i> restricting the val- 
ues ofc so that (c)), //(Ai rz s(c)) < 1/2. These are sets constructed from grid 
boxes whose corresponding entry in the W^ M ' N ^ vectors is above a certain value. 

3. Definei](c) = argmin c , gM \fi(A^(c))—fj,(A^, T Jcf))\. Given a value ofc, 77(c) determines 

the set A^, T z Jr)(c)) that best matches the p-measure of Aj(c), as required by item 3 of 
Definition [7| 

4- Set c* = argmax c6R p M (A+(c), A^,, (77(c))). The value ofc* is selected to maximise the 
coherence. 

5. Define A z := A+(c*) and A^ z) := A^ z) ( V (c*)). 
Remark 6. 

1. One can repeat Algorithm^ replacing Aj(c) and A^,Jc) with A~(c) = Uri4 / ( M - iV )(t)<c ^ 
and AzJc) = \J i:W (M,N)^ TjZ ^ <c Bi respectively See [20] for further details. 

2. Care should be taken regarding the sign of W^ M,N '(z) and W^ M,N '(^(r, z)). Visual 
inspection may be required in order to check that the vectors have the same "parity" . 

7.2 Coherent Sets for a ID discrete time nonautonomous system 

We return to the map cocycle $ and Perron-Frobenius cocycle described in Example [1] and 
identify coherent sets. We use two methods: firstly, inspection of the composition of maps 
as perturbations of maps with invariant sets, and secondly using the general method of 
Algorithm 5. 

The map cocycle <£> is defined in terms of a map H a which has an almost-invariant set, 
and this gives rise to a family of coherent sets in the following way. Recall the definitions 
of the maps T iy i = 1, . . . , 4 and shift space Q determined by the adjacency matrix B. The 
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maps Ti have the property that if Bij = 1, then any inner R factors cancel in Tj o T{. More 
generally, for any u 6 fl, we have cancellation of all intermediate R factors: 

u, -) = R s o # Vfe _ Mo o ■ • • o H auo o R~\ (31) 
where s, t G {0, 1} are given by 

s(u>,k) = \ ?' ^- l0dd ' and t( U ,k) = \ ?' 
v y \ 1, uneven, v y \ 1, w > 2. 

For the map ifo ; the interval [0, 0.5] is invariant. Moreover, [0, 0.5] is almost-invariant for 
H a with p M ([0, 0.5]) = 1 - 2a. By QH), we if we set 

i CTfcw = iT^ fc )([0, 0.5]), for each fceN, (32) 

then 

P/j,(A a k u , A a k+iu) = 1 - 2a Wfc . 

Thus < A CT fe w ^ is a family of p -coherent sets with p = 1 — 2max{ai, . . . ,a 4 } = 0.843. 
I J fceN 

In the same way, the invariant set [0.5,1] of Hq leads to a family |i? s ^' fc ^([0.5, l])} fcgN of 
Po-coherent sets with the same po- 

In order to demonstrate the methods of this article, we now show how Algorithm [5] can 
be used. We may use the Oseledets subspaces computed in Section 15.41 to find a family of 
coherent sets. First we apply Algorithm [5] to find a coherent set for the time step k = to 

k — 1. We calculate p^ \ A+(c), A+^c) j as c varies over the elements of the vector Z^ 20 ' 10 ^); 

see Figure IT21 (left). The maximum value of p M (A+(c), A+ W (r](c))j is 0.890. The set is 
found to be the interval [0.11,0.58] of length ijl{AJ) = 0.47; see Figure IT21 (right). 




-1.0 -0.5 0.0 0.5 1.0 -i , 



Figure 12: (left): The function p M ^4+(c), A+ u (r](c))J takes its maximum on the interval 
(—0.352, —0.176) and so we take the midpoint c* = —0.264 as the optimal threshold, (right): 
Taking this optimal threshold (shown in dashed green) for the eigenvector f^ 0,1 ^ (oj) identifies 
the coherent set = [0.11,0.58] (shown in dark orange). 

We note that the set Aj found by Algorithm [5] is not the same as the A^ produced 
by the intuitive construction (I321) . In the latter case, A w = [0,1/2], A au > = [0,1/2], and 



25 



p^Aaj, A auj ) = 1 — 2a UJo = 1 — 2a\ = 1 — n/20 « 0.843, significantly lower than the value of 
0.890 found using Algorithm |5j 

We may extend Algorithmic in order to find a sequence of coherent sets {A alU j}f =Q . Since 
we require the measure of a sequence of coherent sets to be constant, we seek to maximize 
the mean value of over a given time range as we vary the measure of the sets. 

Algorithm 6 (To determine a sequence of maximally coherent sets over a range of 
times co, ... , a K uj). 

1. Follow steps 1.-3. of Algorithm^ for each k = 0, . . . , K — 1 using r = 1 to obtain sets 
A\ (c). 

2. Let c k {£) := argmin cgK Jc)) - 

3. Compute t := argmax te(0 5] £ (i+ Jc^)), i+ fc+laJ (c fc+ i(f))) . 
^ F O rfc = 0,...,A:-l, define A a M u :=A+. u (c k (e*)). 

To demonstrate Algorithm [6], we use the approximate Oseledets functions f2(cr k U)), k = 
0, . . . , 5, to find a sequence of six coherent sets {Ar^j/Lo for the map cocycle <£>. Plot- 
ting \Y!=oPv (A% u (c k (fy),A* h+lu (c k+1 (£))) a g ainst ^ ( see Figure USD, we find a unique 
maximum of 0.891, which occurs at £* = 0.47. 




Figure 13: The graph of p M := jELo/V (^%uj( c k{^)) , ^k+i^Ck+i^))) against £, where 
we take A+ ku} (c k (£)) for i < 0.5 and (c k (£)) otherwise. The maximum 0.891 occurs at 
£* = 0.47. The red section of the curve corresponds to A\ (c k (£)) and the blue section to 
A- h Jc k (£)). 

Figure [8] shows the graph of f^ ,w \a k u) with the threshold c k (£*) for k = 0, . . . , 5, and 
in each case the set A^ kui (c k (£*)) is indicated by shading. Since coherent sets are required to 
be connected, we must find the interval closest to each A+ k ^(c k (£*)). For k = 0, 2, 3, 4, 5 the 
set A+ kuj (c k {£*)) is itself an interval and we set A a k u = A+ k ^(c k (£*)). The set A+ U (c k (£*)) has 
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k 





1 2 3 


4 5 




1 


2 3 2 


4 3 


A i 


[0.11,0.58] 


[0.12,0.59] [0.35,0.82] [0.07,0.54] 


[0.35,0.82] [0.35,0.82] 


Pk 


0.89 


0.87 0.87 0.96 


0.90 0.87 



Table 1: Coherent sets A^k^ and the values of p M (A^^, A^+i^) for k = 0, . . . , 5. 



two components, [0.12.0.58] and [0.60,0.61], and so we set A auj = [0.12,0.59]. Table [T] lists 
the coherent sets A^k^ and the values of p M (A^k^, A a k+T. u ) for k = 0, . . . , 5. 

It is interesting to compare the locations of the coherent sets with their correspond- 
ing maps in the mapping cocycle, see Figure HH As for the previously constructed family 




T a u] — T2 
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T„2,,, = T-? 
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1 




A 



Xt 5 w — T 3 




Figure 14: Graphs of T a k u showing A a k u x A^+i^ for k — 0, 



5. 



|i? s (^' fc )([0.5, l])} fcgN , the coherent sets alternate between two positions separated by a rota- 
tion of approximately 0.25. However, the mean value of p M is greater for the sequence A a k u 
constructed from Algorithm [6] since in each case the coherent set matches up well with local 
maxima and minima of the preceding map. 



7.3 Coherent Sets in a 2D continuous time nonautonomous system 

We apply Algorithm [5] to the Oseledets subspaces W {m > m) (z) and W/ (80 ' 40) (£(10, z)) cal- 
culated in Section 15.5.11 and displayed in Figure [TTJ The optimal coherent sets A+ and 
£(10 z) are obtained at the threshold values c* = 0.0043 and i](c*) = 0.0052, which gives 
p At (v4+(c*), At 10 ^(r)(c*))) = 0.9605, see Figure [T5l For A~ and i« 10j j the optimal threshold 
is at c* = -0.0040 and r]{c*) = -0.0051 and p^(i"(c*), i" (lM (7/(c*))) = 0.9599. The coher- 
ent sets at Af and Af (1Q , and the images of sample points in Af are shown in Figure dBJ 
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T 




-0.015 -0.01 -0.005 0.005 0.01 0.015 



0^ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ■— 1 

-0.1 -0.05 0.05 0.1 

Figure 15: Thresholding curve p M (4j(c), A~ (lM (77(c))) and (p M (i+(c), i+ 10 2) (r/(c))) are 
plotted in grey and black, respectively. The optimal threshold is marked with a rectan- 
gle. 

In Figure [16] we note that the grey set on the left at time z = (0, 1, 1.5) flows 
approximately to the light grey set A^t 10 z \ on the right at time £(10, z). Similarly for the 
black sets A~ and AZ 1Qz y This carrying of the time z coherent sets to the time £(10, z) 
coherent sets by the aperiodic flow is only approximate, as p^A^c*), A^, lQ Jr)(c*))) = 
0.9605 and p fM (A~(c*), AZ 10 z Jrj(c*))) = 0.9599. Thus, we expect a loss of about 10% under 
the advection of the flow. Figure [16] also zooms onto A^ and A~ to demonstrate this loss of 
mass. To make this loss even more apparent, we continue to flow forward for 50 time units. 
Figure [T7] shows that the (black) coherent sets A+ and A~ do indeed disperse over time, 
however at a much slower rate than the arbitrarily chosen (grey) sets. The coherent sets 
A+, A~ are just single elements of a time parameterised family {A% t AZ t z y$t>o of coherent 
sets that at any given initial time describe those sets that will disperse most slowly over a 
duration of 10 time units. 

8 Final Remarks 

We have formulated a new mathematical and algorithmic approach for identifying and track- 
ing coherent sets in nonautonomous systems. Our new approach generalises existing suc- 
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Figure 16: [Top left] The coherent sets (grey) and A~ (black) and [Top right] A^ 10 2 ) 
(grey) and A« 10 > (black) are identified by thresholding the Oseledets functions. Overlays 
of 0(10, z, A+) (grey) and 0(10, z, A~) (black) are also shown. [Bottom left] Zooms of A + (z) 
and A~(z) . [Bottom right] Overlays of 0(10, z, Af) (grey/black dots) on A^, 1Qz s (white), 
displaying the loss of mass over 10 time units duration from z = (0, 1, 1.5). 
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Figure 17: Trajectories of the perturbed system (1241) for e — 1. The large light grey (A^ 
and black (AZ t z - ) ) blobs are the coherent sets identified by our approach. The other (medium 
grey) blobs are chosen nearby the coherent sets to show strong mixing away from the coherent 
regions. 
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cessful transfer operator methodologies that have been used in the autonomous setting. Our 
constructions address the question raised by [30J of how to study strange eigenmodes and 
persistent patterns observed in forced fluid flows in the general time-dependent situation. 
Future work will include applying these techniques to detect and track mobile coherent re- 
gions in oceanic and atmospheric flows, extending significantly the flow times studied in 
[2U ?] and [?]. 
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